clear, clc, close all;

load level3part01.mat
data = part01;

% load Data20220816154330013all20minlevel3.mat
% data = Data20220816154330013all20minlevel3(:, 1);

% load Data20220812134102013all20minlevel5.mat
% data = Data20220812134102013all20minlevel5;

% begin = 7217000;
% length = 7710660;
% step = 1;
% Y = data(begin: step: begin + length);
Y = data;
X = (1: 1: size(Y, 1))';

plot(X, Y, '-^');

[Yup, Ylow] = envelope(Y, 50, 'rms');

hold on
plot(X, Yup);

% findpeaks(-Yup, 'MinPeakDistance', 400, 'MinPeakProminence',0.2, 'Annotate','extents');
[pks, locs, widths, proms] = findpeaks(Yup, 'MinPeakDistance', 50, 'MinPeakProminence', 0.1);

hold on
plot(locs, pks, '^');

% plot(locs, zeros(size(locs)), '*');
locs_mid = locs(1:end-1) + round(diff(locs) / 2);
locs = locs_mid;
plot(locs_mid, zeros(size(locs)), 's');

% a = locs(8741):10:locs(8742);a = a';b = sin(a);plot(a, b, 'k*');
% figure, plot(Y(locs(199): locs(202)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
segs = cell(size(locs, 1)-1, 1);
segs_data = cell(size(locs, 1)-1, 1);
max_ll = 0; min_ll = 1e8; max_id = 0; min_id = 0;
for i = 1: size(locs, 1) - 1
    segs{i, 1} = Yup(locs(i): 1: locs(i+1));
    segs_data{i, 1} = Y(locs(i): 1: locs(i+1));
    lllll = size(segs{i, 1}, 1);
    if lllll > max_ll
        max_ll = lllll;
        max_id = i;
    end
    if lllll < min_ll
        min_ll = lllll;
        min_id = i;
    end
end

max_lll = 0; min_lll = 1e8; max_idd = 0; min_idd = 0;
for i = 1: size(segs, 1)
    segl = size(segs{i, 1}, 1);
    if segl >= min_ll * 2
        nnnn = floor(segl / min_ll);
        segs{i, 1} = segs{i, 1}(1: nnnn: end);
    end
    llllll = size(segs{i, 1}, 1);
    if llllll > max_lll
        max_lll = llllll;
        max_idd = i;
    end
    if llllll < min_lll
        min_lll = llllll;
        min_idd = i;
    end
end

nn = 2 ^ ceil(log2(max_lll));
ll = nn/2 + 1;
% fft_feats = cell(size(segs, 1), 1);
fft_feats = zeros(size(segs, 1), ll);
for i = 1: size(fft_feats, 1)
    xx = segs{i, 1};
    yf = fft(xx, nn);
    P2 = abs(yf/nn);
    P1 = P2(1: ll);
    P1(2:end-1) = 2*P1(2:end-1);
%     P1 = P2(1: ll) .^ 2;
%     P1 = P2(1: ll) * 2;
    fft_feats(i, :) = P1;
%     figure; plot(P1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[coeffs, ~, latent, ~, explained] = pca(fft_feats);
zz = fft_feats * coeffs(:,1:2);

% [indices, C] = kmeans(zz, 8);
indices = dbscan(zz, 0.002, 5);
figure
plot(indices(indices~=-1))
% plot(indices)

figure; 
gscatter(zz(:, 1), zz(:,2), indices)
axis equal